# 01_NYMBY_OLS.R
# Created: 2021-10-4 Taka-aki Asano
# Last Modified: 2024-02-16

# package
require("readr")
require("dplyr")
require("tidyr")
require("miceadds")
require("lme4")
require("clubSandwich")


# load dataset
data <- read_csv("HLWconjoint20210330.csv")


# processing of data
data2 <- dplyr::select(data, 
                       no, time, rl, hlw_su, hlw_fsu, psm01, psm02, psm03, 
                       psm04, psm05, psm06, psm07, psm08, psm09, psm10, 
                       psm11, psm12, psm13, psm14, psm15, psm16, 
                       sex, age, ac_ba, mar, group, risk_g, q4_1_1)
data2$hlw_fsu <- ifelse(data2$hlw_fsu == data2$rl, 1, 0)


# variables
## PSM
data2$PSM <- rowSums(data2[,paste0("psm", c("01", "02", "03", "04", "05", "06", 
                                            "07", "08", "09", "10", "11", "12", 
                                            "13", "14", "15", "16"))])
summary(data2$PSM)
data2$PSM_Group <- ifelse(data2$PSM >= median(data2$PSM), "HPSM", "LPSM")
data2$PSM_Group <- factor(data2$PSM_Group, levels = c("LPSM", "HPSM"))
aggregate(PSM ~ PSM_Group, data2, mean)

## APS
data2$APS <- rowSums(data2[,paste0("psm", c("01", "03", "12", "15"))])
summary(data2$APS)
data2$APS_Group <- ifelse(data2$APS >= median(data2$APS), "HAPS", "LAPS")
data2$APS_Group <- factor(data2$APS_Group, levels = c("LAPS", "HAPS"))
aggregate(APS ~ APS_Group, data2, mean)

## CPV
data2$CPV <- rowSums(data2[,paste0("psm", c("04", "06", "07", "08"))])
summary(data2$CPV)
data2$CPV_Group <- ifelse(data2$CPV >= median(data2$CPV), "HCPV", "LCPV")
data2$CPV_Group <- factor(data2$CPV_Group, levels = c("LCPV", "HCPV"))
aggregate(CPV ~ CPV_Group, data2, mean)

## COM
data2$COM <- rowSums(data2[,paste0("psm", c("02", "09", "10", "16"))])
summary(data2$COM)
data2$COM_Group <- ifelse(data2$COM >= median(data2$COM), "HCOM", "LCOM")
data2$COM_Group <- factor(data2$COM_Group, levels = c("LCOM", "HCOM"))
aggregate(COM ~ COM_Group, data2, mean)

## SS
data2$SS <- rowSums(data2[,paste0("psm", c("05", "11", "13", "14"))])
summary(data2$SS)
data2$SS_Group <- ifelse(data2$SS >= median(data2$SS), "HSS", "LSS")
data2$SS_Group <- factor(data2$SS_Group, levels = c("LSS", "HSS"))
aggregate(SS ~ SS_Group, data2, mean)

## Female
data2$Female <- ifelse(data2$sex == 2, 1, 0)

## Undergraduate degree
data2$Undergraduate <- ifelse(data2$ac_ba == 5, 1, 0)

## Graduate degree
data2$Graduate <- ifelse(data2$ac_ba == 6, 1, 0)

## Married
data2$Married <- ifelse(data2$mar %in% 2:4, 1, 0)

## Having children
data2$Children <- ifelse(data2$mar %in% 3:4, 1, 0)

## Rural
data2$Rural <- ifelse(data2$group %in% 4:6, 1, 0)


# regression
## Model1
reg1 <- lm.cluster(
  hlw_su ~ PSM + Female + age + Undergraduate + Graduate + 
    Married + Children + risk_g + q4_1_1 + Rural, 
  data = data2, cluster = "no"
)
summary(reg1)

## Model2
reg2 <- lmer(
  hlw_su ~ PSM + Female + age + Undergraduate + Graduate + 
    Married + Children + risk_g + q4_1_1 + Rural + (1 | no),
  data = data2
)
summary(reg2)
coef_test(reg2, vcov = "CR2", test = "Satterthwaite")

## Model3
reg3 <- lm.cluster(
  hlw_su ~ APS + CPV + COM + SS + Female + age + Undergraduate + 
    Graduate + Married + Children + risk_g + q4_1_1 + Rural, 
  data = data2, cluster = "no"
)
summary(reg3)

## Model4
reg4 <- lmer(
  hlw_su ~ APS + CPV + COM + SS + Female + age + Undergraduate + 
    Graduate + Married + Children + risk_g + q4_1_1 + Rural + (1 | no),
  data = data2
)
summary(reg4)
coef_test(reg4, vcov = "CR2", test = "Satterthwaite")
